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00 i ABSTRACT 



Context. Twilight studies have proved to be important tools to analyze the atmospheric structure with interesting consequences on the 

characterization of astronomical sites. Active discussions on this topic have been recently restarted in connection with the evaluation of Dome 

C, Antarctica as a potential astronomical site and several site-testing experiments, including twilight brightness measurements, are being 

prepared. 

Aims. The present work provides for the first time absolute photometric measurements of twilight sky brightness for ESO-Paranal (Chile), 

which are meant both as a contribution to the site monitoring and as reference values in the analysis of other sites, including Dome C. 

Methods. The UBVRI twilight sky brightness was estimated on more than 2000 FORSI archival images, which include both flats and standard 

Q stars observations taken in twilight, covering a Sun zenith distance range 94°-l 12°. 

5-H Results. The comparison with a low altitude site shows that Paranal V twilight sky brightness is about 30% lower, implying that some fraction 

C/j , of multiple scattering has to take place at an altitude of a few km above the sea level. 
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1 . Introduction community for what seems to be the new frontier of ground- 
based astronomy, i.e. Dome C - Antarctica. This site is excep- 

The quality of an astronomical site is determined by several pa- ^-^^^^ ^^^^^ ^^^^ ^^^p^^^^ 3^^^^^^ ^^ extremely good seeing 

rameters, which may vary according to the wavelength range of conditions reported by Lawrence et al. (-^OOTt, several stud- 
interest. For the optical and near-IR domain, these include typ- -^^ j^^^^ ^^^^^^ ^^^^ 1^^ ^^^^^^^ ^f precipitable water vapor 
ical seeing, sky transparency, number of clear nights, humid- ^j^^^j^^ ^^^pj^^ ^^ ^ j^^ ^y.^ emission, could imply that this is 
ity, night sky brightness, amount of precipitable water vapor, ^^e best site for IR and sub-millimetric ground-based astron- 
dust and aerosols. While the seeing, extinction, sky brightness ^^^ ^ ^^^-^^ ^^^^^ ^1^^ characteristics of Dome C has been 
and other quantities are commonly measured at most observa- ^g^enfly presented by Kenyon & Storey C005}. 
tories, the twilight brightness is not. This is mainly because 

the relevant information on the typical atmospheric conditions One of the main concerns is related to the high latitude 

can be derived from other measurements obtained during the of this site. This, in fact, causes a significant reduction in 

night. Nevertheless, twilight observations provide an indepen- the amount of dark time with respect to equatorial observa- 

dent tool for probing the overhead atmosphere under much tories, thus posing some doubts about the effective exploita- 

higherfluxconditions, thus allowing more accurate results. The tion of the exceptional seeing in the optical. The possibility 

interested reader can find an extensive review on this topic in of opening spectral windows otherwise unaccessible from the 

the classical textbook by Rozenberg (1966 1. ground is in itself a valid and sufficient scientific driver for 

Very recently, the twilight has received particular attention Dome C. Nevertheless, arguments in favor of Dome C as a site 

due to the growing interest of the international astronomical for optical astronomy have been put forward. Among these, a 

smaller average contribution by the scattered moonlight to the 

Send ojfprint requests to: F. Patat global background and a cleaner atmosphere have been advo- 

* Based on observations made with ESO Telescopes at Paranal cated as features that may possibly compensate for the reduced 

Observatory. dark time (Kenyon & Storey 2005) . In particular, since the last 
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phases of twilight (also known as deep twilight) are dominated 
by multiple scattering (Rozenberg 1 1 966> . the amount of scat- 
tered sunlight is strongly dependent on the amount of aerosol 
in the lower atmospheric layers (see for example Ougolnikov et 
al.2004i. In a supposedly low aerosol content site like Dome 
C, this effect is expected to be very low and this, in turn, would 
allow one to start the observations earlier than at normal sites. 
Even though the argument has a good physical ground, di- 
rect on site measurements are still lacking. In this respect, it 
is worth mentioning that a couple of dedicated experiments 
for sky brightness measurements are currently being setup (A. 
Moore, J. Storey 2006, private communications). 

In spite of the large number of investigations done in the 
past in this field, absolute twilight brightness measurements 
are rather rare, especially for large observatories placed in top 
rated sites. To our knowledge, the only published work on twi- 
light observations in the Johnson-Cousins standard system is 
the one by Tyson & Gal ("1993^ who however, given their pur- 
poses, report only uncalibrated data for CTIO. In the light of 
these fact, both with the purpose of characterizing Paranal also 
from this new point of view and to provide the community with 
absolute reference values obtained over a large time baseline, 
we present here for the first time UBVRI twilight sky bright- 
ness measurements. 

The paper is organized as follows. In Sec.|2we introduce 
the basic concepts through a simplified model (which is dis- 
cussed in more detail in Appendix A), while in Sec. |3] we 
describe the observations, data reduction and calibration. The 
UBVRI twilight sky brightness measured at ESO-Paranal is 
presented and discussed in Sec. |3 while Sec. |5] summarizes 
the results obtained in this work. 

2. The twilight problem 

The calculation of scattered flux during twilight is a rather 
complicated problem that requires a detailed treatment of mul- 
tiple scattering (see for example Blattner et al. 119741 and 
Wu & Lu 1988 1 and an accurate description of the atmo- 
spheric composition and the physical phenomena taking place 
in the various layers (Divari & Flo tnikova 119661 Rozenberg 
[1966 1. Notwithstanding the large amount of work done in the 
'60 and in the '70, the problem is still matter of investiga- 
tions (see for example Anderson & Lloyd 1990 Ougolnikov 
[T9991 Ougolnikov & Maslov 2002 Ekstrom 2002 Ougolnikov, 
Postylyakov & Maslov 2004, Postylyakov 2004, Mateshvili et 
al. 12)05 1. While it is well beyond the purposes of the present 
work to explore the problem from a theoretical point of view, 
we deem it is interesting to introduce a simple single-scattering 
model, which is useful both to understand the basic principles 
of twilight and to provide a quick comparison to the observed 
data. The assumptions and the model itself are discussed in 
Appendix|X] to which we refer the interested reader for the de- 
tails, while here we concentrate on the model predictions only. 
The calculated zenith {a-Q) UBVRI sky brightness as a 
function of Sun zenith distance ^, computed using the average 
broad band extinction coefficients for Paranal (Patat 2003a i, 
are plotted in Fig. ^ As one can immediately see, the single 
scattering component drops below the night sky brightness be- 
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Fig. 1. Model twilight sky brightness at zenith. The thick curves 
include the night sky contribution, while the thin lines indicate 
the scattered component only. The vertical dashed-dotted line 
marks the Sun zenith distance when the lower boundary layer 
height is 120 km. The upper scale indicates the lower Earth's 
boundary layer height in km. 



tween ^=99° and ( =100°, indicating that from this point on 
multiple scattering is the only contributor to the observed flux, 
as shown by Ugolnikov & Maslov i2002 ) on the basis of polar- 
ization measurements. 

Another aspect to be considered is the fact that the tran- 
sition to the flatter part of the atmospheric density profile 
(see Fig. lA.2> definitely occurs during the multiple scattering- 
dominated phase. Since multiple scattering takes place with 
higher probability where the density is higher, i.e. in the lower 
atmospheric layers, the explanation given by Tyson & Gal 
( 1993 ) for the observed brightness decline rate during twilight 
does not seem to be coiTect. In fact, these authors interpret the 
observed values as the pure consequence of the lower shadow 
boundary height change, neglecting extinction and multiple 
scattering. They conclude that, since their observations have 
been taken when 100< h~ < 400 km (where h, is the height 
of the lower Earth's shadow boundary along the zenith direc- 
tion; see also Appendix A), the sky brightness rate is directly 
related to the slope of the density law in that region of the at- 
mosphere. Nevertheless, the calculation of Sun's ephemeris for 
the site and epoch of Tyson & Gal's observations shows that, 
in the most extreme case (see their Table 1, R filter), it was 1° 
. 8 < ^ < 7?4 (where ip - ( - 90). As the reader can figure out 
from Table IaTI this implies that h, <80 km in all cases, i.e. 
well within the steep part of the density profile. Therefore, the 
fact that the observed rate and the rate expected from pure sin- 
gle scattering in the higher atmospheric layers are consistent, 
is just a coincidence. 
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Fig. 2. Distribution of twilight observations in Alt-Az coordi- 
nates for TSF (filled symbols) and LSE (empty symbols). The 
astronomical azimuth has been replaced with the difference in 
azimuth between the telescope pointing and the Sun, Aa©. 

3. Observations, Data Reduction and Calibration 

In order to measure the twilight sky brightness on Paranal, we 
have used archival calibration data obtained with the FOcal 
Reducer/low dispersion Spectrograph (hereafter FORSl), 
mounted at the Cassegrain focus of ESO-Antu/Melipal 8.2m 
telescopes (Szeifert 2002). The instrument is equipped with a 
2048x2048 pixels (px) TK2048EB4-1 backside thinned CCD 
and has two remotely exchangeable collimators, which give 
a projected scale of 0'.'2 and O'.'l per pixel (24yum x 24yum). 
According to the used collimator, the sky area covered by the 
detector is 6.'8x6'8 and 3'4x3.'4, respectively. For this study we 
have selected only the data obtained with the lower resolution 
collimator and the 4-port high-gain read-out mode, since this 
combination is the most used for imaging with FORS 1 . With 
this setup the read-out noise is 5.5 electrons (e"). 

For our purposes, we have selected two sets of data. The 
first is composed by broad-band UBVRI twilight sky flats 
(hereafter TSF), which are regularly obtained as part of the cal- 
ibration plan. In the current implementation, after taking a test 
exposure the observing software estimates with a simple algo- 
rithm the integration time for the first exposure in a series of 
4 frames. Subsequent exposure times are adjusted on the ba- 
sis of the previous exposure level and this allows one to ob- 
tain high signal-to-noise images with a rather constant counts 
level, which is typically around 20,000 ADUs. This is achieved 
with exposure times which range from 0.25 sec up to 5 min- 
utes. Given these values, the sensitivity of FORSl in the vari- 
ous passbands and the typical twilight sky brightness behaviour 
(see for example Fig.^, these observations are expected to ap- 
proximately cover the Sun zenith distance range 94° < ^ < 
101°, i.e. still within the nautical twilight. Since the most im- 
portant part of this analysis concerns the deep twilight, it is 



Fig. 3. Zenith twilight sky surface brightness in the U pass- 
band from TSF (empty symbols) and LSE (filled symbols) with 
\a\ <40°. The vertical dashed lines mark the end of nautical 
(left) and astronomical (right) twilight. The solid curve traces 
the simple model described in the text. The horizontal lines 
indicate the U dark time night sky brightness average value 
(dashed) for Paranal (Patat .2003aj and the +lcr levels (dotted). 



clear that an additional set of data is required to complement 
the sky flats. 

The calibration plan of FORSl includes the observation 
of standard stars fields (Landolt"1992"i in UBVRI passbands, 
which are regularly taken during twilight, typically just after 
the sky flats sequence is completed. For calibration purposes, 
a fraction of these exposures are obtained using relatively long 
integration times (typically 40 sec for U and 20 sec for BVRI) 
which, at an 8 m-class telescope, are suflicient to bring the sky 
background at exposure levels which are suitable to our pur- 
poses. In fact, the bulk of these observations covers the range 
100° < ^ < 112°, i.e. well into astronomical twilight. For the 
sake of clarity we will indicate them as long exposure standards 
(LES). 

In order to collect a statisticaUy significant sample, we have 
retrieved from the ESO Archive all suitable TSF obtained in 
UBVRI passbands from 01-01-2005 to 30-09-2005, for a to- 
tal of 1083 frames {U: 148, B: 208, V: 226, R: 261, /: 240). 
Since a much higher night-to-night spread is expected in the 
deep twilight phase due to the natural fluctuations of the night 
sky emission and also because LSE are less frequently obtained 
than TSF, a larger time interval must be considered. To this aim, 
we have put together the LES sample collecting all suitable im- 
ages obtained from 01-01-1999 (i.e. shortly after the beginning 
of FORSl operations) to 30-09-2005, thus covering almost 6 
years, for a total of 3388 frames {U: 923, B: 611, V: 609, R: 
635, /: 610). All images were processed within the xccdred 
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package of IRAF'. Due to the large amount of data and the 
purpose of this work, the bias subtraction was performed using 
a pre-scan correction only, while flat-fielding was achieved us- 
ing a stack of all TSF in each given passband as a master flat, 
which was then adopted to correct each TSF and LES frame. 

Due to the nature of the data, there is no need for taking 
care of the possible presence of crowded stellar fields or bright 
extended objects, as it is the case for night sky brightness es- 
timates (see Patat '2003bV For this reason, the background in 
each image was measured using a simple and robust mode 
estimator. To avoid possible vignetting and flat fielding prob- 
lems, only the central 1024x1024 pixels were considered. On 
these spatial scales and due to improper flat-fielding, FORS 1 is 
known to show variations of the order of a few percent, while 
the gradients in the twilight sky are much smaller (see Chromey 
& Hasselbacher 1996 1 and can be safely neglected. Therefore, 
the mode < / > of pixel intensity distribution is assumed as the 
best estimate of the sky background. For each filter, this is con- 
verted into a surface brightness in the Johnson-Cousins system 
via the foUowing relation 

b - -2.5 log < / > -1-2.5 \og{p^texp) + niQ 

where p is the pixel scale (arcsec px"'), texp is the exposure 
time (seconds) and mo is the instrumental photometric zero- 
point for the given passband. No color correction has been ap- 
plied, for two reasons: a) colour terms in FORS 1 are very small 
(Patat l200'3ab : b) color correction depends on the intrinsic color 
which, for the twilight sky, changes according to time and posi- 
tion. Given this and the fact that night-to-night fluctuations are 
the dominating source of noise in the measurements, the color 
correction can be safely neglected. As for the instrumental ze- 
ropoints, we have used the average values computed for the 
period of interest using data taken during photometric nights 
only. Variations in the zeropoints are known to take place (Patat 
I2003at . mainly due to the aging of telescope reflective surfaces 
but, again, these are smaller than the inherent sky variations. 
Finally, no airmass correction has been applied and this is ex- 
pected to produce an additional spread in the data close to the 
end of astronomical twilight and a systematic increase in the 
average sky brightness in that region. Due to the large number 
of pixels used (A^ >10^), the typical internal photometric error 
is expected to be less than 1%. 

4. Twilight sity brightness at ESO-Paranal 

Not being obtained specifically for twilight brightness mea- 
surements, the data are inhomogeneously distributed on the 
sky. In fact, in an Alt-Az plot where the ordinary azimuth is 
replaced by the difference in azimuth between the sky patch 
and the Sun (Aoq = a - Gq), the data points tend to cluster in 
two regions, which correspond to evening and morning obser- 
vations (Fig. |2j. Besides target azimuth and altitude, for each 
data point we have computed a series of other quantities, which 



' IRAF is distributed by the National Optical Astronomy 
Observatories, which are operated by the Association of Universities 
for Research in Astronomy, under contract with the National Science 
Foundation. 
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Fig. 4. Same as Fig.|3]for the B passband. 

Table 1. Twilight sky brightness fitted parameters in the range 
95° < ^ < 105°. All values are expressed in mag arcsec"^. 



Filter 



flo 



deg"' 



02 



rr-2 



y 

deg"' 



u 


11.78 


1.376 


B 


11.84 


1.411 


V 


11.84 


1.518 


R 


11.40 


1.567 


I 


10.93 


1.470 



-0.039 0.24 1.23±0.01 

-0.041 0.12 1.24±0.01 

-0.057 0.18 1.14±0.02 

-0.064 0.29 1.09+0.03 

-0.062 0.40 0.94±0.03 



are relevant for the subsequent analysis. These include Sun az- 
imuth and altitude. Sun-target angular separation. Moon phase. 
Moon altitude and Moon-target angular separation. To avoid 
contamination from scattered moonlight in the LSE, we have 
selected only those data points for which the Moon was below 
the horizon. 

Since the twilight sky brightness for a given Sun zenith dis- 
tance changes with the position on the sky, in order to study its 
behavior as a function of ( it is necessary to make a selection 
on the Alt-Az coordinates. Given the nature of the available 
data, which appear to be rather concentrated (Fig. IS, it seems 
reasonable to restrict the analysis to zenith region only. In or- 
der to have a sufficient amount of measurements, we have used 
all data points with zenith distance \a\ <40°, which is of course 
expected to cause some additional spread in the observed re- 
lation. The results are presented in Figs. 13171 for the UBVRI 
passbands respectively. 

As expected, the single scatter model drops much faster 
than the actual observations. Whilst for VRI the model devi- 
ates from the data around ( -97°, for B and especially for 
U the model underestimates the surface brightness already at 
( <96°, indicating that multiple scattering is more efficient at 
shorter wavelengths. This is in agreement with the findings of 
Ougolnikov & Maslov C2002j . who have shown that the con- 
tribution of single scattering in the phases immediately follow- 
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Fig. 5. Same as Fig.0]for the V passband. 



Fig. 7. Same as Fig.|3]for the / passband. 
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Fig. 6. Same as Fig.|3for the R pass band. 



ing sunset is about 40%, 60%, 70% and 80% in U, B, V and 
R respectively. These fractions remain roughly constant un- 
til ( <95°, after which the role of single scattering becomes 
weaker and weaker and multiple scattering takes rapidly over 
In all passbands, the night sky brightness level is reached at 
around f=105°-106°. 

In order to give a more quantitative description of the obser- 
vations, we have fitted the surface brightness data in the range 
95° < ( < 105° using second order polynomials of the form 
ao + ai(^ -95) + a2(^ -95)^ , with ^ expressed in degrees and the 
surface brightness in mag arcsec"^. The results are presented in 
Table [0 where we have reported also the RMS deviation from 
the fitted function cr and the slope y deduced from a linear fit to 
the data in the range 95° < ^ < 100°, i.e. during the interval typ- 



Fig. 8. Broad band zenith twilight sky colors. The curves have 
been computed using second oder polynomials fitted to the 
observed data. For comparison, the colors of the Sun are 
U - B=0.13, B - y=0.65, V - R^O.52 and V - 7=0.81, while 
those of the night sky at Paranal are t/ - B=-0.36, B-V -1.03, 
V - R^0J4 and V - 7=1.90 (Patat l2003a> . 



ically used to obtain TSF exposures, when the contribution by 
the night sky is still moderate. A first aspect to be noticed is the 
spread shown by the data points around the mean laws, which 
are due to the night-to-night variations in the atmospheric con- 
ditions. The dispersion becomes particularly large in the 7 pass- 
band, where the fluctuations appear to be quite pronounced. As 
for the decay rate during nautical twilight, we notice that this 
tends to decrease for increasing values of wavelength. 
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To convert the values of y reported in Table [Qinto surface 
brightness variation per unit time, one has to multiply them by 
difldt, which is given by: 

d(fi dHg 

— = cos (fi sin Hq cos Sq cos 4> 

dt dt 

where Hq and 6q are the hour angle and declination of the 
Sun respectively, while (p is the site latitude (see for example 
Smart'W??!. Since dH^ldt ^0.25 deg min"', for tp ~0 (i.e. at 
the time of sunrise and sunset) for Paranal (0=-24?6) one ob- 
tains d(p/dt-Q.23 deg min"' at the equinoxes (6q=Q). Applying 
this factor to the values reported in the last column of Table[n 
one obtains brightness decline rates that range from 0.28 (U) 
to 0.22 (/) mag arcsec"^ min '. These values can be directly 
compared with those reported by Tyson & Gal ('1993V The 
data available to these authors did not allow them to quantify 
the differences between the various filters and hence they re- 
port a mean value of y=0.23+0.02 mag arcsec"^ min"' which 
is within the range defined by our data. 

With the aid of the second order best fit relations, we have 
computed the color curves presented in Fig.|8l Due to the dis- 
persion of the observed data, these colors have to be regarded 
only as indicative, especially in the region ^ >102°, where the 
inherent night sky brightness fluctuations start to be signifi- 
cant, particularly in the red passbands. It is interesting to note 
that while U - B and V - R colors do not change very much 
as the Sun sinks below the horizon, significant changes take 
place m B - V and especially in V - /. In principle one ex- 
pects that since multiple scattering boosts the light at shorter 
wavelengths with respect to the pure single scattering compo- 
nent, the overall color gets bluer and bluer as the Sun deep- 
ens below the horizon. Then, at some point, the night sky glow 
which has completely different colours, starts to contribute and 
the colors progressively turn to those typical of the night sky. 
As a matter of fact, the observed U - B, B - V and V - R 
curves indeed show this behaviour (see Fig. |S|l, while V - I 
turns steadily redwards. This is due too the interplay between 
input Sun spectrum, scattering efficiency, extinction, multiple 
scattering and the emergence of the night glow, which combine 
together in quite a complicated way. This is clearly illustrated 
in Fig.|5] where we present twilight spectra obtained on Paranal 
with FORSl and his twin instrument F0RS2. The sky spectra 
were extracted from spectrophotometric standard stars observa- 
tions taken during twilight (see Table |2j and were wavelength 
and flux calibrated with standard procedures in IRAF The ex- 
posure times ranged from 10 sec to 120 seconds and the signal- 
to-noise ratio was increased meshing all pixels in the direction 
perpendicular to the dispersion, after removing the region of 
the detector affected by the well exposed standard star spec- 
trum. For comparison. Fig. |9l shows also the typical dark time 
spectrum of Paranal (Patat 2 003 aJ and the solar spectrum. 

With the exception of the first spectrum, which was ob- 
tained with a very low resolution (~130 A FWHM), the re- 
maining data allow one to detect quite a number of details. For 
if <12°, i.e. during the nautical twilight, the spectrum is rather 
different from that of the Sun, even though it shows clear solar 
features, like the Call H&K lines and the G-band at about 4300 
A. The Rayleigh scattered Sun flux gives a clear contribution 



Table 2. Basic data for twilight sky spectra shown in Fig.|9l 



Date 


UT stai-t 


AA/A 


Alt 


Az 


V 






A(FWHM) 


deg 


deg 


deg 


2001-09-21 


10:06:58 


130 


82.6 


117.5 


6.5 


2003-11-29 


09:02:25 


13 


65.8 


122.9 


9.4 


2003-02-28 


09:45:52 


16 


27.5 


299.3 


11.9 


1999-04-16 


09:54:36 


12 


41.1 


2334 


14.8 



in the region bluewards of 5000A down to (p-l5°, after which 
the pseudo-continuum of the night sky emission takes over 

For ifi <9°, the contribution by night sky emission 
lines is very weak. Characteristic lines like [OI]5577A and 
[OI]6300,6364 A (the latter overimposed on a O2 absorption 
band) are barely visible, while the OH Meinel bands start 
to appear above 8000 A. One remarkable exception are the 
Nal D lines, which are known to be present during the so 
called sodium flash (Rozenberg 1 1 966> . A similar phenomenon 
is present for the [OI]6300,6364 A doublet (see for example 
Roach & Gordon [l973> . which is indeed visible in Fig. |9l al- 
ready at ip-9°4. With the onset of astronomical twilight, the 
spectrum in the red is more and more dominated by the OH 
bands. Another remarkable fact is the behaviour of the H2O 
(6800, 7200 A) and O2 (7600 A) molecular absorption bands. 
During the bright twilight, when single scattering is still rel- 
evant {ip=6°5), they akeady appear to be significant, but they 
become even deeper in the multiple scattering-dominated phase 
(^=9?4), due to the longer optical path traveled by the multiply 
scattered photons. For higher values of ip they progressively 
disappear due to the weakening of the scattered Sun's contin- 
uum. 

In order to estimate the effects of site altitude on the twilight 
sky brightness, we have compared the results presented here 
with data obtained at the Southern Laboratory of the Sternberg 
Astronomical Institute (Moscow, Russia) during three morning 
twilights on December 9-11, 2002 (see Fig.llO>. This facility is 
located within the Crimean Astrophysical Observatory (CrAO), 
at a latitude of 44?7 North and 600 m a.s.l.. The observations 
were performed using a wide-field CCD-camera with a field 
of view of 8°x6° and exposure times that ranged from 0.01 to 
18 seconds. Photometric calibration was achieved using field 
stars included in the Tycho-2 Catalog (Hog et al. I2000> . The 
photometric passband of this instrument is fairly similar to the 
Johnson-Cousins V, with a colour correction of the order of 
0.01 mag for the (B - V) colour range shown by twilight data. 

The two data sets clearly show that the twilight back- 
ground at CrAO is systematically brighter; the difference is 
constant during the dark twilight period and vanishes at night- 
fall. More precisely, the comparison between V band data 
in moonless conditions at CrAO (zenith) and ESO-Paranal 
(\a\ < 20°) shows that the mean difference in the Sun de- 
pression range 5°5 < ip < 11?0 is AV = 0.27+0.03. On the 
other hand, the typical atmospheric pressure value for ESO- 
Paranal is P 1-743 hPa, and for CrAO during the observations 
P2=961 hPa. Interestingly, the magnitude difference implied 
by the pressure ratio at the two sites is -2.5 log(Pi/P2)=0.28, 
which is indeed very similar to the measured difference AV. 
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Fig. 9. Twilight spectra obtained at Paranal at different Sun de- 
pression angles. For convenience, the spectra have been cor- 
rected to the corresponding V flux shown by the photometric 
data. Spectral resolution and sky patch diff'er from spectra to 
spectra (see Table |2j- For presentation the spectra have been 
vertically shifted by the following amounts: -1-0.50 {ip-9.A), 
+035 {ip^n.9), -0.10 ((^=14.8), -1.0 {tp >=18). The dotted 
line traces the solar spectrum, normalized to the continuum of 
the ip-9A spectrum at 4500 A. 

Therefore, we can conclude that the deep twilight sky bright- 
ness is proportional to the atmospheric pressure or, equiva- 
lently, to the atmospheric column density above the observer 
In turn, this implies that the light undergoes multiple scatter- 
ing throughout the whole atmosphere and not only in the upper 
layers. Given that the difference in altitude between Paranal and 
CrAO is only 2 km, the observations we present here suggest 
that some fraction of multiple scattering has to take place in the 
first few km above the sea level. 

5. Discussion and Conclusions 

In this paper we have presented for the first time absolute 
UBVRI twilight brightness measurements for the ESO-Paranal 
Observatory (Chile) spanning almost 6 years. These measure- 
ments will serve as reference values for the similar studies 
which will be soon conducted at Dome C, Antarctica, as part 
of the site testing campaign. The planned in situ spectropho- 
tometric measurements will finally clarify whether this excep- 
tional location shows a lower aerosols content, as expected both 
due to the icy soil and to its large distance from the sea coast 
(Kenyon & Storey 2005 1. 

The twilight sky brightness measurements presented here 
were obtained from VLT-FORS 1 archival data not specifically 
taken for this kind of studies. Also, due to the large telescope 
diameter, the initial twilight phases (0° < ip <6°) are not cov- 
ered. In a sense this is quite unfortunate, since for these low Sun 
depression angles the lower shadow's boundary passes through 



Fig. 10. Comparison between the zenith twilight sky bright- 
ness measured on Paranal (small symbols) and Crimean 
Astrophysical Observatory (large symbols) for the V passband. 



the atmospheric layers below ~30 km (see Table lA.ll . where 
the ozone and aerosol stratospheric concentration is maximum. 
These phases are in fact used to retrieve ozone and aerosol pro- 
files, using both intensity and polarization measurements (see 
for example Wu & Lu ,1988. Ugolnikov et al. 2004; Mateshvili 
et al. 2005 1. Nevertheless, during the deep twilight, when the 
direct Sun radiation illuminates only the upper atmospheric 
layers and single scattering on air molecules becomes pro- 
gressively less important, the amount of aerosols and ozone 
plays a relevant role through multiple scattering. Therefore, 
even though of much more complicated interpretation, deep 
twilight observations may still give some insights on the con- 
ditions in the lower atmosphere. Moreover, in the context of 
the discussion about the supposedly shorter twilight duration 
at Dome C (see Kenyon & Storey 2005 ), what really matters 
is the behavior during the deep twilight. An example of this 
kind of analysis is shown in Fig. ^2 where the data obtained 
at Paranal are compared to the MCC++ model calculations for 
a site at 2.6 km a.s.l. (see Postylyakov 2004 for a detailed de- 
scription). This code treats the radiative transfer in a spherical 
atmosphere including Ozone, Aerosol and molecular scatter- 
ing, also taking into account the backscattering by the Earth's 
surface. As one can see the overall behavior is fairly well repro- 
duced. The deviations are possibly due to the differences in real 
and model aerosol, since multiple scattering is very sensitive to 
it. The model adopts a urban microphysical model for the first 
10 km of the atmosphere and this is certainly different from 
what is expected for a desertic area close to the sea, as is the 
case of Paranal. Dedicated instruments for twilight sky bright- 
ness monitoring coupled to detailed modeling may indeed give 
a useful contribution to the already existing site testing tools, 
providing independent indications about the overhead aerosol 
profile. 
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Fig. 11. Comparison between Paranal zenith sky brightness 
and the MCC++ model (Postylyakov 20Q4j. The dashed curve 
traces the pure model solution without the contribution of night 
sky emission. 

Some interesting results are obtained comparing the es- 
timates obtained for Paranal (2600 m a.s.l.) with those of a 
significantly lower site like CrAO (600 m a.s.l.). Despite the 
fact that the bright twilight and night sky brightnesses are very 
close at the two sites, during the deep twilight Paranal is about 
30% darker than CrAO in the V passband (see Fig. I10> . Due 
to the higher altitude, Paranal suffers from a lower extinction 
which, if all other atmospheric properties are identical and mul- 
tiple scattering takes place mostly in the troposphere (5-10 km, 
Ougolnikov & Maslov 2002 1, would then turn into a brighter 
twilight sky. The observations actually show the opposite be- 
haviour and the brightness ratio is fairly consistent with the 
atmospheric pressure ratio (see previous section). A natural in- 
terpretation of this fact is that a fraction of the multiple scat- 
tering events takes place at heights which are lower than it was 
originally thought, say below 3 km from the sea level. 

Whether this is due to the lower density of air molecules, to 
a smaller amount of ground level aerosols or to a combination 
of the two needs further investigation and the comparison with 
other astronomical sites at even higher altitudes, like Mauna 
Kea. 

Appendix A: A simple semi-analytical model for 
twilight brightness 

A.1. Basic assumptions 

The model is based on the following simplifying assumptions: 
(i) Earth is a sphere with radius /?o=6380 km; (ii) the atmo- 
sphere extends up to A/?=400 km and the numerical density 
n{h) of the scattering particle density is given by the MSIS-E- 
90 model profile (Hedin .l991j ; (iii) the effect of atmospheric 
refraction can be neglected; (iv) the Sun is a point source and 



all incoming sun rays are parallel; (v) only single scattering 
plus attenuation is considered; (vi) Rayleigh scattering by air 
molecules is the only source of sunlight diffusion. While some 
of these assumptions are reasonable, (v) and (vi) are a bit crude 
and will certainly lead to discrepancies between the model re- 
sult and the actual observations. For a more sophisticated single 
scattering model taking into account refraction and the pres- 
ence of ozone and aerosols, the reader is referred to Adams, 
Plass & Kattawar J 1 974( 1 . 

A.2. Model geometry 

As we have said, we will assume that the Earth is a sphere of ra- 
dius Rq and that the observer is placed at an elevation hs above 
the sea level. Since we will consider only small values of hs (<3 
km), we make the further simplifying assumption that the hori- 
zon is a plane tangent to the sphere of radius Rq + hs in O, (see 
Fig. lA.U . i.e. neglecting the horizon depression. We will indi- 
cate with (fi the Sun depression (ip > 0) and with a the zenith 
distance of the generic sky patch under consideration. For the 
sake of simplicity we will derive the sky brightness only along 
the great circle passing through the zenith (Z) and the Sun. This 
angular distance is counted positively in the direction of the 
Sun, so that negative angles indicate sky patches in the anti- 
Sun direction. Under these simplifying assumptions, the lower 
boundary of the Earth's shadow is described by a straight line, 
which is tangent to the sphere in Hq, and which is indicated by 
a dotted-dashed line in Fig. IA.il When the observer is looking 
into the generic direction a, the corresponding line of sight will 
cross the lower shadow boundary in Pq and the contribution to 
the observed flux will come from all the scattering elements 
along the segment PqPi, where Pi indicates the intersection 
between the line of sight and the upper atmospheric boundary, 
which is placed at an altitude AR above the sea level^. 

In this geometry, a volume element placed in P receives the 
sun light, which is attenuated along its path OP, and it scatters 
the photons into the observer's direction, with a scattering an- 
gle = n/2 - a + if. according to the scattering phase function 
0(0), which obeys to the usual normalization condition 



f 

J4n 



(D(6») d£l = 1 



(A.l) 



Before reaching the observer, it undergoes the extinction 
along OP. The total flux will finally result from the integration 
along the illuminated segment P^P]. In order to compute the 
required quantities, we must first derive a series of useful rela- 
tions. Since one has to evaluate the particle density along the 
generic light path, it is fundamental to know the height h - GP 
above the sea level for any given point along the trajectory. 
This is particularly simple for the travel between O and Pi, for 
which one can easily derive the following relation: 



hop^ - ^jP- + 2l(Ro + hs) cos a + (Pq + /z.,)^ - Rq 



(A.2) 



where / is the coordinate along OPi (l-O in O). Another 
useful relation is the one that gives the lower limit for the in- 
tegral along the line of sight, i.e. the optical path /q between O 

- The upper limit is set just for numerical reasons. 
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Fig. A.l. Geometry of the problem. For the sake of clarity, the 
observer's elevation hs has been exaggerated. 



and P(). In order to get the required expression, one first needs 
to find the length of HqPq and ho, i.e. the height above the sea 
level of Pq. These can be obtained after some trigonometry and 
are respectively 



HqPq — Rq 


tani^- 


sin a /I - cos ip 
cos(a - ip)\ cos if 


Ro 


and 







/io= yJHoPl+Rl-Ro 

Having these two equations at hand, one can easily find Iq: 



/o = ^(Ro + h,)^ cos2 a + hl-hl + 2Ro{h - hs) 
-{Ro + hs) cos a 

As for the upper limit h along the line of sight, this turns 
out to be: 

h = ^J(Ro + hsf- cos2 a + IRaiAR - hs) + AR^ - h^^ 
-{Rq + hs) cos a 

The next step is the calculation of Hqp, the height along the 
sun ray. In order to do so, we introduce the perigee height 6, i.e. 
the minimum distance between the Earth's surface and the sun 
ray passing through P, which can be expressed as a function of 
the / coordinate along PqPi as follows: 

6 = (I - Iq) cos{a - ip) 

Having this in mind, one can derive the length of QH as 
follows: 

QH^ y/lRoiAR - 6) + AR^ - 6^ 



Fig.A.2. Normalized density (solid line) and temperature 
(dashed line) profiles according to the MSIS-E-90 model 
(Hedinn"991^. For comparison, the dots trace the values of the 
US Standard Atmosphere (McCartnev .1976. Tab. 2.6). 

If we then introduce a coordinate q along QP (with q=Q in 
Q), considering the fact that (QH - qf + (Ro + 6)^ - (Rq + hf, 
we finally obtain 



V= ^{QH-qf + {Ra + 6)^-Ro 



(A.3) 



with < ^ < QP. The optical path of the unscattered 
sun ray QP is given by the sum of QH and HP, the latter be- 
ing HP - H\iP() - (I - Iq) sin(a - tp). By means of Equations 
IA.2I and IA.3I one can now compute the height above the sea 
level (and hence the particle number density) along the light 
path. The height along the zenith direction, h^ - TZo, can be 
readily derived and it is given by 



h, = Ro 



1 - cos (p 
cosip 



From this one can figure out how fast h, grows when the 
sun depression (p increases (see also Tab. lA. 11 1. This, coupled to 
the rapid decrease of the atmospheric density as a function of 
height, is the cause for the very quick drop in the twilight sky 
surface brightness. 

A.3. Density profile and opticai depth 

As we have anticipated, for the density profile we have adopted 
the MSIS-E-90 model^ density profile (Hedin [79911 1. which is 
presented in Fig. IA.2I As the plot shows, the global profile 
can be roughly described by two laws: one exponential (for 
h <120 km, the so called homosphere) and a power law (for 
120 < h < 400 km, the so called thermosphere). Clearly, as the 
Sun depression increases, the lower Earth's shadow boundary 

^ http://modelweb.gsfc.nasa.gov/models/msis.html 
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Table A.l. Height of lower boundary of Earth's shadow at 
zenith (h^) and at 60° from zenith along the solar meridian in 
the Sun's direction. 



by the observer per unit time, unit area and unit wavelength is 
given by 



V 


hz 


ho{a = 60) 


(deg) 


(km) 


(km) 


3 


8.8 


8.0 


6 


35.1 


29.9 


9 


79.5 


63.3 


12 


142.5 


106.7 


15 


225.1 


159.1 


18 


328.3 


220.1 



will pass through more and more tenuous atmospheric layers 
and, therefore, the change in density slope should turn into a 
change in the twilight sky surface brightness decline. Under the 
assumptions made for this simplified model, this should happen 
when if ^101°, i.e. close to the end of nautical twilight. 

Since in the following section we will be interested in the 
product between the number density and the extinction cross 
section, we can derive this, for a given passband, assuming the 
measured extinction coefficient k{A) (in mag airmass"') and in- 
tegrating the previous density profile along the vertical (i.e. at 
airmass 1). In fact, assuming that all the extinction comes from 
Rayleigh scattering, one can write: 

T,(/l) = no Cf vrW) I dh 

Jh, no 

Then, considering that t,(A)=I.Q^6k(A), the product 
no Cevf(/i) can be readily derived from the previous relation 
and used in all further optical depth calculations. 

A.4. Scattering cross section and ptiase function 

The /I '* wavelength dependency of the Rayleigh scattering 
cross section is implicitly taken into account by the extinc- 
tion coefficients, which are to be considered as input data to 
the model and not as free parameters. As for the scattering 
phase function we have used the canonical expression for air 
molecules (McCartnev I 1 976> : 



<D(6i) = [1 +0.9324 cos2(0)] 

lOTT 

where ffie multiplicative constant comes from the normal- 
ization condition expressed by Eauation lA.il 



A. 5. Scattered flux 

The scattered flux can be computed in the same way as done, 
for example, by Krisciunas & Schaefer ( 1991 1 for the Moon 
light. If Lq is the luminosity of the sun at a given wavelength 
(expressed in photons per unit time and per unit wavelength), 
ffie flux received by the Earth at the top of the atmosphere is 
Fg - Lq/AttcI^, where d-l AU. If we now consider an in- 
finitesimal volume element dV placed along the line of sight at 
a distance / from O, ffie number of scattered photons received 



df. 



0^-T(QP) 



n^ 



n[h(l)] CJA) 



P 



(OP) 



dV 



(A.4) 



Given the geometry of the problem, the infinitesimal vol- 
ume element can be written as dV - dS dl = nf'^^ dl, where 
(f) is the semi-amplitude of the angle subtended by dS at the 
distance of the observer. Since the solid angle subtended by dS 
is dQ. - jKJp', the surface brightness produced by ffie volume 
element is simply 

db^^= Fj^e-"(2P) „[/j(/)] c^„(A) <D(0) e-'^'^''> dl 

and the total surface brightness is finally obtained integrat- 
ing along the line of sight within the illuminated region: 






dbdl 



Finally, to take into account the contribution by the night 
sky emission, we have added to the computed flux the one im- 
plied by the average values measured for Paranal (/zj=2.6 km) 
and reported by Patat (i2003a» Tab. 4). With this, the model is 
completely constrained and ffiere are no free parameters. 
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